

El problema del churn es un problema que afecta a todas las compañías, pero en especial a las Telcos, con una tasa superior al 30% debido a la alta competencia entre ellas y la gran facilidad de cambiar de una a otra compañía. Según un estudio de Daemon Quest, retener a un cliente cuesta entre cinco y quince veces menos que captar a uno nuevo. Es por ello por lo que las telcos, al igual que otras compañías, están poniendo mucho interés en la retención de los clientes.
El objetivo de este NoteBook es predecir el comportamiento de los clientes para su posterior retencion en la compañia. Para ello se analizarán todos los datos relevantes de los clientes descargados en el siguiete enlace:

Esta descripción del data set es dada por el enunciado, posteriormente analizaremos paso a paso los datos del enunciado.
Cada fila representa un cliente, cada columna contiene atributos de cliente descritos en la columna Metadata. Los datos brutos contienen 7043 filas (clientes) y 21 columnas(variables/atributos). La Columna “Churn” es nuestro target.

Cargamos el dataset y librerias para su posterior manipulación:
import pandas as pd
import numpy as np
churn =pd.read_csv("WA_Fn-UseC_-Telco-Customer-Churn.csv")
churn.info()
Como podemos observar coinciden tanto el numero de columnas como el numero de filas que tenemos en el dataset. Pero el tipo columna no esta definido como tal.
Posteriormente procedemos a pasar las columnas correspondientes de tipo object a tipo float. Pero antes de ello vamos a explorar un poco el contenido de cada columna y si estas contienen valores nulos.
# Contenido de cada columna
churn.apply(set,axis = 0)
# Conteo de valores nulos
churn.isnull().sum()
Porcentaje de custumer-churn para saber si es posible con el dataset generar la prediccion o debereiamos hacer un tratamiento para generar mas custumer-churn
churn['Churn'].value_counts(sort=True,normalize = True)

Tras un pequeño analisis, coincidimos con el enunciado del dataset. El siguiente paso es diferenciar nuestras varibales numericas de las categoricas, y realizar los cambios de tipologia necesarios:
col_id = ['customerID']
col_churn = ["Churn"]
col_numericas = ['tenure','MonthlyCharges','TotalCharges']
col_categoricas = ['gender', 'SeniorCitizen', 'Partner', 'Dependents','PhoneService', 'MultipleLines',
'InternetService','OnlineSecurity', 'OnlineBackup', 'DeviceProtection', 'TechSupport',
'StreamingTV', 'StreamingMovies', 'Contract', 'PaperlessBilling','PaymentMethod']

Tras intentar realizar la siguiente operación:
churn[var_numericas] = churn[var_numericas].astype(float)
Nos muestra un error, el cual nos indica que dentro de estos numeros hay una string. Si hacemos la vista atras, podemos obvservar que la unica variable que no estaba en formato float era TotalCharges.
(churn.loc[:,['customerID','TotalCharges']]).sort_values(by =['TotalCharges']).head()
Para encontrar estos valores hemos ordenado los valores, obligandole a colocar en primer lugar los valores de caracter string y por ultimo los tipo float
churn[churn.TotalCharges == ' '].loc[:,['TotalCharges']].count()
Estos valores respresentan un 0.15% del dataSet, con lo cual procedemos a eliminarlos.
churn = churn[churn.TotalCharges != ' ']
churn = churn.reset_index()[churn.columns]
churn[col_numericas] = churn[col_numericas].astype(float)
Una vez eleminado estos outliers, empecemos un analisis gráfico con boxplots y la distribución de los datos
import matplotlib.pyplot as plt
import seaborn as sns
fig, ax = plt.subplots()
fig.set_size_inches(5, 7)
ax = sns.boxplot(x="Churn", y="tenure", data=churn,width=0.3, notch=True)
Si hubieramos hecho un analisis sin la diferencia de churn-custumer no hubieramos visto estos outliers dentro de yes-churn. Como estos datos son escasos no podemos permitirnos el lujo de quitarnoslos, con lo cual vamos a suavizarlos poniendo como toque maximo, el punto maximo de nuestro percentil. Hacemos esto debido a que son muy pocos y ademas no se salen demaisado del marco.
churn['tenure'] = np.where((churn.Churn == 'Yes')&(churn.tenure > 67),69,churn.tenure)
fig, ax = plt.subplots()
fig.set_size_inches(5, 7)
ax = sns.boxplot(x="Churn", y="tenure", data=churn,width=0.3, notch=True)
fig, ax = plt.subplots()
fig.set_size_inches(5, 7)
ax = sns.boxplot(x="Churn", y="MonthlyCharges", data=churn,width=0.3, notch=True)
fig, ax = plt.subplots()
fig.set_size_inches(5, 7)
ax = sns.boxplot(x="Churn", y="TotalCharges", data=churn,width=0.3, notch=True)
churn[(churn.Churn == 'Yes')&(churn.TotalCharges > 6800)]
Como hemos dicho anteriormente eliminar estos datos no es lo correcto debido a que escasea informacion para los casos de custumer - churn es positivo. Lo que vamamos hacer es suavizar estos datos para que no entorpecer el aprendizaje. Para este caso vamos a colocarle la media de los valores que se salen de lo normal.
churn[(churn.Churn == 'Yes')&(churn.TotalCharges > 6800)].loc[:,['TotalCharges']].mean()
churn['TotalCharges'] = np.where((churn.Churn == 'Yes')&(churn.TotalCharges > 6800),7337.694792,churn.TotalCharges)
fig, ax = plt.subplots()
fig.set_size_inches(5, 7)
ax = sns.boxplot(x="Churn", y="TotalCharges", data=churn,width=0.3, notch=True)
Hemos descubierto que hay un pequeño salto entre los outliers, vamos a utilizar ese salto para cortar.
churn['TotalCharges'] = np.where((churn.Churn == 'Yes')&(churn.TotalCharges > 6500),6500,churn.TotalCharges)
import seaborn as sns
sns.pairplot(churn[col_numericas + col_churn ], kind="scatter", diag_kind= 'kde',height=4, hue="Churn")
# Cambiamos una varible binaria a categórica para poder hacer el plot
churn['SeniorCitizen'] = churn['SeniorCitizen'].replace({1:"Yes",0:"No"})
for i in col_categoricas:
sns.pairplot(churn[col_numericas + col_categoricas ], kind="scatter", diag_kind= 'kde',height=4, hue=i)
Una de las conclusiones a destacar tras el plot, es la modificación de las variables categoricas: ['OnlineSecurity','OnlineBackup','DeviceProtection','TechSupport','StreamingTV', 'StreamingMovies']
Que tienen como variables 'Yes', 'No', 'No internet services', las cuales podemos reducirlas a 'Yes', 'No'. Hemos decidido hacer esto debido a la cantidad y distribucíon de los datos.

Procederemos a la exploración de las variables categoricas, tomando como primera acción el cambio mencionado anteriormente y convirtiendo las varibales a numéricas:
churn[col_categoricas] = churn[col_categoricas].replace({"No internet service":"No"})
churn[col_categoricas] = churn[col_categoricas].replace({"No phone service":"No"})
churn[col_categoricas] = churn[col_categoricas].replace({"Female":1,"Male":0})
churn[col_categoricas] = churn[col_categoricas].replace({"Yes":1,"No":0})
churn[col_churn] = churn[col_churn].replace({"Yes":1,"No":0})
# Realizamos un Get dummies para las variables con mas de dos categorias
churn_dummies = pd.get_dummies(churn[col_numericas + col_categoricas + col_churn])
churn = (churn.loc[:,['customerID']]).join(churn_dummies)
col_categoricas = churn.nunique()[churn.nunique() == 2].keys()
col_categoricas
import plotly.offline as py
import plotly.graph_objs as go
import plotly.tools as tls
import plotly.figure_factory as ff
py.init_notebook_mode(connected=True)
dat_rad = churn_dummies[col_categoricas]
def plot_radar(df,aggregate,title) :
data_frame = df[df["Churn"] == aggregate]
data_frame_x = data_frame[col_categoricas].sum().reset_index()
data_frame_x.columns = ["feature","yes"]
data_frame_x["no"] = data_frame.shape[0] - data_frame_x["yes"]
data_frame_x = data_frame_x[data_frame_x["feature"] != "Churn"]
#contamos los Yes que son 1
trace1 = go.Scatterpolar(r = data_frame_x["yes"].values.tolist(),
theta = data_frame_x["feature"].tolist(),
fill = "toself",name = "Yes",
mode = "markers+lines",
marker = dict(size = 5)
)
#contamos los No que son 0
trace2 = go.Scatterpolar(r = data_frame_x["no"].values.tolist(),
theta = data_frame_x["feature"].tolist(),
fill = "toself",name = "No",
mode = "markers+lines",
marker = dict(size = 5)
)
layout = go.Layout(dict(polar = dict(radialaxis = dict(visible = True,
side = "counterclockwise",
showline = True,
linewidth = 2,
tickwidth = 2,
gridcolor = "white",
gridwidth = 2),
angularaxis = dict(tickfont = dict(size = 10),
layer = "below traces"
),
bgcolor = "rgb(243,243,243)",
),
paper_bgcolor = "rgb(243,243,243)",
title = title,height = 700))
data = [trace2,trace1]
fig = go.Figure(data=data,layout=layout)
py.iplot(fig)
#Ploteamos
plot_radar(dat_rad,1,"Churn")
plot_radar(dat_rad,0,"Non Churn")

En este apartado generaremos nuevas dimensiones de las variables continuas:
from sklearn.preprocessing import PolynomialFeatures
pf = PolynomialFeatures(degree=2, interaction_only=False,
include_bias=False)
features = pf.fit_transform(churn[col_numericas])
features
Una vez generada una matriz de caracteristicas, miremos con que grados a jugado en cada columna para crear el dataframe
pd.DataFrame(pf.powers_, columns=['tenure_degree','MonthlyCharges_degree','TotalCharges_degree'])
col_features = ['tenure','MonthlyCharges','TotalCharges','tenure_2','tenure*MonthlyCharges','tenure*TotalCharges',
'MonthlyCharges_2','MonthlyCharges*TotalCharges','TotalCharges*2']
churn_features = pd.DataFrame(features, columns = col_features)
# Añadimos algunas medias que nos parecen interesantes
churn_features['MonthlyCharges_mean'] = churn_features['MonthlyCharges'].mean()
churn_features['tenure_mean'] = churn_features['tenure'].mean()
# Genreamos los dos cuerpos de matrices para el posterior modelado
churn_B = (churn.loc[:,['customerID']]).join((churn_features.join(churn[col_categoricas])))
churn_A = churn
En esta matriz hemos dejado las varibles continuas sin manipular, adjuntando sus variables categoricas que anteriormemte hemos tratato.
churn_A.head()
churn_A.info()
En esta matriz hemos jugado con las variables continuas añadiendo medias y combinaciones lineas de estas. Además de las variables categóricas.
churn_B.head()
churn_B.info()

Uno de los fallos mas habituales es dar a la maquina datos no escalados. Esto entorpeceria el aprendizaje de la maquina debido al peso de las varibales de alto valor, dejando aun lado la importancia casi de las variables categóricas. Por ello es necesario realziar un escalado a corde al dataset. Otro punto importantes es ver las correlaciones de las variables para ver si estamos dando la misma información.
from sklearn.preprocessing import StandardScaler
std = StandardScaler()
norm_churn_A = std.fit_transform(churn_A[col_numericas])
norm_churn_A = pd.DataFrame(norm_churn_A,columns=col_numericas)
sns.pairplot(norm_churn_A, kind="scatter", diag_kind= 'kde',height=4)
churn_A_final = (churn_A.loc[:,['customerID']]).join((norm_churn_A.join(churn_A[col_categoricas])))
churn_A_final.info()
#Preparamos el frame de correlaciones
correlation = churn_A_final.corr()
#Extraemos el nombre de las columnas
matrix_cols = correlation.columns.tolist()
#convertimos el frame en array
corr_array = np.array(correlation)
#Visualizacion de datos con plotly
trace = go.Heatmap(z = corr_array,
x = matrix_cols,
y = matrix_cols,
colorscale = [[0.0, 'rgb(165,0,38)'], [0.1111111111111111, 'rgb(215,48,39)'], [0.2222222222222222, 'rgb(244,109,67)'],
[0.3333333333333333, 'rgb(253,174,97)'], [0.4444444444444444, 'rgb(254,224,144)'], [0.5555555555555556, 'rgb(224,243,248)'],
[0.6666666666666666, 'rgb(171,217,233)'],[0.7777777777777778, 'rgb(116,173,209)'], [0.8888888888888888, 'rgb(69,117,180)'],
[1.0, 'rgb(49,54,149)']],
colorbar = dict(title = "Pearson Correlation coefficient",
titleside = "right"
) ,
)
layout = go.Layout(dict(title = "Matriz A de correlaciones",
autosize = False,
height = 720,
width = 800,
margin = dict(r = 0 ,l = 210,
t = 25,b = 210,
),
yaxis = dict(tickfont = dict(size = 9)),
xaxis = dict(tickfont = dict(size = 9))
)
)
data = [trace]
fig = go.Figure(data=data,layout=layout)
py.iplot(fig)
from sklearn.preprocessing import StandardScaler
std = StandardScaler()
norm_churn_B = std.fit_transform(churn_B[col_features])
norm_churn_B = pd.DataFrame(norm_churn_B,columns=col_features)
sns.pairplot(norm_churn_B, kind="scatter", diag_kind= 'kde',height=4)
churn_B_final = (churn_B.loc[:,['customerID']]).join((norm_churn_B.join(churn_B[col_categoricas])))
churn_B_final.info()
#Preparamos el frame de correlaciones
correlation = churn_B_final.corr()
#Extraemos el nombre de las columnas
matrix_cols = correlation.columns.tolist()
#convertimos el frame en array
corr_array = np.array(correlation)
#Visualizacion de datos con plotly
trace = go.Heatmap(z = corr_array,
x = matrix_cols,
y = matrix_cols,
colorscale = [[0.0, 'rgb(165,0,38)'], [0.1111111111111111, 'rgb(215,48,39)'], [0.2222222222222222, 'rgb(244,109,67)'],
[0.3333333333333333, 'rgb(253,174,97)'], [0.4444444444444444, 'rgb(254,224,144)'], [0.5555555555555556, 'rgb(224,243,248)'],
[0.6666666666666666, 'rgb(171,217,233)'],[0.7777777777777778, 'rgb(116,173,209)'], [0.8888888888888888, 'rgb(69,117,180)'],
[1.0, 'rgb(49,54,149)']],
colorbar = dict(title = "Pearson Correlation coefficient",
titleside = "right"
) ,
)
layout = go.Layout(dict(title = "Matriz B de correlaciones",
autosize = False,
height = 720,
width = 800,
margin = dict(r = 0 ,l = 210,
t = 25,b = 210,
),
yaxis = dict(tickfont = dict(size = 9)),
xaxis = dict(tickfont = dict(size = 9))
)
)
data = [trace]
fig = go.Figure(data=data,layout=layout)
py.iplot(fig)
Podemos observar como en la matriz de correlaciones aparece nuestro feature engineering, marcandonos con un color mas azul estas variables correlacionadas.

Antes de empezar a lanzar los modelos debemos separar nuestros datos de entrenamiento y testing.
from sklearn.model_selection import train_test_split
train_A,test_A = train_test_split(churn_A_final, test_size = .25 ,random_state = 111)
train_B,test_B = train_test_split(churn_B_final, test_size = .25 ,random_state = 111)
A continuacion separamos las variables a predecir y col_id que no interviene en nuestro entrenamiento
# Matriz A
cols_A = [i for i in churn_A_final.columns if i not in col_id + col_churn]
train_XA = train_A[cols_A]
train_YA = train_A[col_churn]
test_XA = test_A[cols_A]
test_YA = test_A[col_churn]
# Matriz B
cols_B = [i for i in churn_B_final.columns if i not in col_id + col_churn]
train_XB = train_B[cols_B]
train_YB = train_B[col_churn]
test_XB = test_B[cols_B]
test_YB = test_B[col_churn]
Preparamos una funcion que nos enseñe las carecteristicas importantes de cada entrenamiento, para posteriromente comparar y sacar conclusiones
from sklearn.metrics import confusion_matrix,accuracy_score
from sklearn.metrics import roc_auc_score,roc_curve,scorer
import statsmodels.api as sm
def churn_prediction(algoritmo, training_x, testing_x, training_y, testing_y):
algoritmo.fit(training_x,training_y)
prediccion = algoritmo.predict(testing_x)
probabilidad = algoritmo.predict_proba(testing_x)
conf_matrix = confusion_matrix(testing_y,prediccion)
model_roc_auc = roc_auc_score(testing_y,prediccion)
fpr,tpr,thresholds = roc_curve(testing_y,probabilidad[:,1])
print (algoritmo)
print ("Accuracy Score : ",accuracy_score(testing_y,prediccion))
print ("Area bajo la curva : ",model_roc_auc,"\n")
#Curva ROC
grafico1 = go.Scatter(x = fpr,y = tpr,
name = "ROC : " + str(model_roc_auc),
line = dict(color = ('rgb(22, 96, 167)'),width = 2),
)
grafico1B = go.Scatter(x = [0,1],y=[0,1],
line = dict(color = ('rgb(205, 12, 24)'),width = 2,
dash = 'dot'))
#plot confusion matrix
grafico2 = go.Heatmap(z = conf_matrix ,x = ["Not churn","Churn"],
y = ["Not churn","Churn"],
showscale = False,colorscale = [[0.0, 'rgb(165,0,38)'], [0.1111111111111111, 'rgb(215,48,39)'], [0.2222222222222222, 'rgb(244,109,67)'],
[0.3333333333333333, 'rgb(253,174,97)'], [0.4444444444444444, 'rgb(254,224,144)'], [0.5555555555555556, 'rgb(224,243,248)'],
[0.6666666666666666, 'rgb(171,217,233)'],[0.7777777777777778, 'rgb(116,173,209)'], [0.8888888888888888, 'rgb(69,117,180)'],
[1.0, 'rgb(49,54,149)']],name = "matrix",
xaxis = "x2",yaxis = "y2"
)
layout = go.Layout(dict(title="Caracteristicas del modelo" ,
autosize = False,height = 600,width = 900,
showlegend = False,
plot_bgcolor = "rgb(243,243,243)",
paper_bgcolor = "rgb(243,243,243)",
xaxis = dict(title = "Ratio falso positivo",
gridcolor = 'rgb(255, 255, 255)',
domain=[0, 0.6],
ticklen=5,gridwidth=2),
yaxis = dict(title = "Ratio verdadero positivo",
gridcolor = 'rgb(255, 255, 255)',
zerolinewidth=1,
ticklen=5,gridwidth=2),
margin = dict(b=200),
xaxis2=dict(domain=[0.7, 1],tickangle = 90,
gridcolor = 'rgb(255, 255, 255)'),
yaxis2=dict(anchor='x2',gridcolor = 'rgb(255, 255, 255)')
)
)
data = [grafico1,grafico1B,grafico2]
fig = go.Figure(data=data,layout=layout)
py.iplot(fig)

Para este apartado se a decidido utilizar los modelos mas utilizados hoy en dia, los cuales son:
from sklearn.linear_model import LogisticRegression
LogisticRegression = LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
intercept_scaling=1, max_iter=100, multi_class='ovr', n_jobs=1,
penalty='l2', random_state=None, solver='liblinear', tol=0.0001,
verbose=0, warm_start=False)
churn_prediction(LogisticRegression, train_XA, test_XA, train_YA, test_YA)
churn_prediction(LogisticRegression, train_XB, test_XB, train_YB, test_YB)
from sklearn.ensemble import RandomForestClassifier
RandomForestClassifier =RandomForestClassifier(bootstrap=True, class_weight=None, criterion='entropy',
max_depth=3, max_features='auto', max_leaf_nodes=None,
min_impurity_decrease=0.0, min_impurity_split=None,
min_samples_leaf=1, min_samples_split=2,
min_weight_fraction_leaf=0.0, n_estimators=100, n_jobs=1,
oob_score=False, random_state=None, verbose=0,
warm_start=False)
churn_prediction(RandomForestClassifier, train_XA, test_XA, train_YA, test_YA)
churn_prediction(RandomForestClassifier, train_XB, test_XB, train_YB, test_YB)
from sklearn.naive_bayes import GaussianNB
GaussianNB = GaussianNB(priors=None)
churn_prediction(GaussianNB, train_XA, test_XA, train_YA, test_YA)
churn_prediction(GaussianNB, train_XB, test_XB, train_YB, test_YB)
from xgboost import XGBClassifier
XGBClassifier = XGBClassifier(base_score=0.5, booster='gbtree', colsample_bylevel=1,
colsample_bytree=1, gamma=0, learning_rate=0.9, max_delta_step=0,
max_depth = 7, min_child_weight=1, missing=None, n_estimators=100,
n_jobs=1, nthread=None, objective='binary:logistic', random_state=0,
reg_alpha=0, reg_lambda=1, scale_pos_weight=1, seed=None,
silent=True, subsample=1)
churn_prediction(XGBClassifier, train_XA, test_XA, train_YA, test_YA)
churn_prediction(XGBClassifier, train_XB, test_XB, train_YB, test_YB)

Tras visualizar los graficos de los modelos, llegamos a la conclusion que el mejor algoritmo para este caso es la LogisticRegression, para la matriz B, dandonos los siguientes resultados:
Accuracy Score : 0.8100113765642776 Area bajo la curva : 0.7311836090903239
Cabe descatar que el feature engeniering ha mejorado los resultados respecto a la matriz A, muy poco, pero la mejora ha sido buena. Recalco, que este resultado puede mejorarse si se estudia la posibilidad de mayores combinaciones lineales entre las variables continuas. Además, los algoritmos que hemos probado, han sido con parametros sencillos, podría caber la posibilidad que los algoritmos como el XGBoost o el RamdomForest superen este score jugando con sus hyperparametros.